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Abstract: In this paper, we derive a closed form solution for recovering the motion of 
an observer relative to a planar surface directly from image brightness derivatives. We 
do not compute the optical flow as an intermediate step, only the spatial and temporal 
intensity gradients at a minimum of 8 points. We solve a linear matrix equation for the 
elements of a 3x3 matrix. The eigenvalue decomposition of its symmetric part is then 
used to compute the motion parameters and the plane orientation. 
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The problem of determining rigid body motion and surface structure from image data 
has been the topic of many research papers in the area of machine vision [1-22]. Many 
approaches based on, tracking feature points [5,11,19,20] or contours [9], using motion 
flow field [1,3,4,10,12,16,17,21,22], texture [2], or image intensity gradients [14,15] have 
been proposed in the literature. 

In the feature point matching schemes, information about a finite number of well- 
seperated points is used to recover the motion (general 8-point 2-frame algorithms of 
Longuet-Higgins [11], Tsai and Huang [20], Buxton et al. [5], and the algorithm of Tsai, 
Huang and Zhu [19] for planar surfaces). These methods require identifying and matching 
feature points in a sequence of images. The minimum number of points required depends 
on the number of image frames. With 2 frames, in most cases, a minimum of 5 points 
results in a unique solution from a set of nonlinear equations. However, using 8 points, 
as in algorithms cited above, one only solves linear equations. Here, it is assumed that 
the more difficult problem of establishing point correspondence has already been solved. 
In general, this involves determining corners along contours using iterative searches. For 
images of smooth objects, it is difficult to find good features or corners. 

For the general case of smooth surfaces, Longuet-Higgins and Prazdny [11] suggested 
a method that uses the optical flow and its first and second derivatives at a single point. 
Later, Waxman and Ullman [21] developed this into an algorithm for recovering the struc- 
ture and motion parameters from a set of nonlinear equations. Subbarao and Waxman 
[17] recently found a closed form solution to the original formulation in [21] for planar 
surfaces. These methods while mathematically elegant are very sensitive to errors in the 
optical flow data since second order derivatives of noisy data are used. 

At the expense of more computation, more robust algorithms have been suggested 
using the optical flow at every image point [1,3,4]. 

Longuet-Higgins [12] has presented a closed form solution for planar surfaces, very 
similar to ours, using the coefficients of the second order optical flow equations. However, 
it is assumed that both components of the flow field have already been computed for a 
minimum of 5 image points. 

By representing a planar surface in the form of a closed contour, Kanatani [9] has 
shown that the surface and motion parameters can be computed by measuring "diam- 
eters" of the contour using line and surface integrals. Here, no point correspondence is 
required. 

Assuming that the planar surface has a uniform texture density, Aloimonos and Chou 
[2] have presented a procedure for computing the motion and surface orientation from 
texture. 

In much of the research work in recovering surface structure and motion from the 
optical flow field, it is assumed that a reasonable estimate of the full optical flow field 
is available. In general, the computation of the local flow field exploits a constraint 
equation between the local intensity changes and the two components of the optical flow. 
However, this only gives the component of the flow in the direction of the intensity 
gradient. To compute the full flow field, one needs additional constraints such as the 
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heuristic assumption that the flow field is locally smooth [7,8]. This, in many cases, leads 
to optical flow fields that are not consistant with the true motion field. 

In an earlier paper, we presented an iterative scheme for recovering the motion of an 
observer relative to a planar surface directly from the image brightness derivatives, and 
without the need to compute the local flow field [14,15]. Further, using a compact vector 
notation, we showed that, at most, two interpretations are possible for planar surfaces 
and derived the relationship between them. Here, we present a closed form solution to the 
same problem. We first solve a linear matrix equation for the elements of a 3x3 matrix 
using intensity derivatives at a minimum of 8 non-colinear points. The special structure 
of this matrix allows us to compute the motion and structure parameters very easily. 



2. Preliminaries 

We first recall some details about perspective projection, the motion field, the brightness 
change constraint equation, rigid body motion and planar surfaces. This we do using 
vector notation in order to keep the resulting equations as compact as possible. 

2.1. Perspective Projection 

Let the center of projection be at the origin of a Cartesian coordinate system. Without 
loss of generality we assume that the effective focal length is unity. The image is formed 
on the plane 2=1, parallel to the xy-plane, that is, the optical axis lies along the z-axis. 
Let R be a point in the scene. Its projection in the image is r, where 

1 A 



Rz 

The 2-component of r is clearly equal to one, that is r • z = 1. 

2.2. Motion Field and Optical Flow 

The motion field is the vector field induced in the image plane by the relative motion of 
the observer with respect to the environment. The optical flow is the apparent motion of 
brightness patterns. Under favourable circumstances the optical flow is identical to the 
motion field (moving shadows or uniform objects in motion could create descrepencies 
between the motion field and the optical flow. Here, we assume that the motion flow field 
and the optical flow are the same) . The velocity of the image r of a point R is given by 

dr d 1 

T> 

—r- — ~~rzz — r-"" 

dt dtK-z 

For convenience, we introduce the notation r^ and R* for the time derivatives of r and 
R, respectively. We then have 

rf= R^ Rf -(R^ (Rf -* )R ' 
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which can also be written in the compact form 

(It • Z) 

since a x (b x c) = (c • a)b — (a • b)c. The vector r* lies in the image plane, and so 
(rf • z) = 0. Further, r^ = 0, if R f || R, as expected. 
Finally, noting that R = (R • z)r, we get 

rf-r-r(zx(R,xr)). 



2.3. Rigid Body Motion 

In the case of the observer moving relative to a rigid environment with translational 
velocity t and rotational velocity w, we find that the motion of a point in the environment 
relative to the observer is given by 

R( = -uxR-t. 

Since R = (R ■ z)r, we can write this as 

R t = -(R- z)wxr-t. 

Substituting for R* in the formula derived above for r^, we obtain 

r t = -(z x (r x (r x u - ^-^rt))) . 
v R- z ' 

It is important to remember that there is an inherent ambiguity here, since the same 
motion field results when distance and the translational velocity are multiplied by an 
arbitrary constant. This can be seen easily from the above equation since the same 
image plane velocity is obtained if one multiplies both R and t by some constant. 

2.4. Brightness Change Equation 

The brightness of the image of a particular patch of a surface depends on many factors. 
It may for example vary with the orientation of the patch. In many cases, however, it 
remains at least approximately constant as the surface moves in the environment. If we 
assume that the image brightness of a patch remains constant, we have 

dE 

-d7 = °> 

or 

dE dr dE _ 

~&7 ' U + ~dt ~ °' 

where dE/dr = (dE/dx,dE/dy,0) T is the image brightness gradient. It is convenient 
to use the notation E r for this quantity and Et for the time derivative of the brightness. 
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Then we can write the brightness change equation in the simple form 

E t ■ r t + E t = 0. 
Substituting for r* we get 

E t -E r -(zx(rx(rxuj- — — t))) = 0. 

XV. * Ml 

Now 

E r -{zx(rx t)) = (E r x z) • (r x t) = {(E r x z) x r) • t, 

and by similar reasoning 

E t • (z x (r x (r x w))) = (((E r x z) x r) x r) • w, 

so we have 

£« - (((#r x z) x r) x r) • w + ^-^r((E r x z) x r) • t = 0. 

To make this constraint equation more compact, let us define c = Et, s = (E r X z) x r, 
and v = — s x r; then, finally, 

C + V • U) -\ ;-S • t = 0. 

R- z 

This is the brightness change equation in the case of rigid body motion. 

2.5. Planar Surface 

A particularly impoverished scene is one consisting of a single planar surface. The equa- 
tion for such a surface is 

Rn = 1, 

where n/ |n| is a unit normal to the plane, and 1/ |n| is the perpendicular distance of the 
plane from the origin. Since R = (R • z)r, we can write this as 

1 
r n = 



Rz' 

so the constraint equation becomes 

c + v • u + (r • n) (s • t) =0. 

This is the brightness change equation for a planar surface. Note again the inherent 
ambiguity in the constraint equation. It is satisfied equally well by two planes with the 
same orientation but at different distances provided that the translational velocities are 
in the same proportions. 
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2.6. Essential Parameters for Planar Surfaces 
The brightness change equation can be written as: 

c + (r x s) • u + (r • n) (s • t) = 0. 
Now, using the identity (r x s) • w = (s x w) • r, we get: 

c + (s x u) • r + (r • n)(s • t) = 0. 

Let us define: 

/ -u>s +w 2 \ 

W = i +w 3 -ui i , 

V — W2 +wi / 



c + r T Ws + r r nt T s = 0, 



then, (s x w) = Ws, in which case, we arrive at 

and after collecting terms: 

c + r T (W + nt T )s = 0. 

Now let: 

f Pi Pi P3 ^ 

P = i P4 Ps P& i = W + nt T , 

V Pi p% P9 J 

We will refer to p t - , i = 1, 2, ..., 9 as the essential parameters (in agreement with Tsai and 
Huang [20]) since these parameters contain all the information about the planar surface 
orientation and motion parameters. Substituting P into the brightness change equation, 
we get: 

c + r T Ps = 0. 

This is a linear constraint equation in terms of the elements of P, and can be used to 
solve for these parameters. We will show how the special structure of P can be exploited 
to recover the motion and structure parameters very easily. 

3. Recovering Essential Parameters 

The nine essential parameters satisfy the following constraint equation: 

c + r T Ps = 0. 

or in terms of p, a vector of length 9 whose i-th element is p»: 

c + a r p = 0, 

where: 

a = (nsi, ns 2 , ris 3 , r 2 si, r 2 s 2 , r 2 5 3 , r 3 si, r s s 2 , r 3 5 3 ) T . 

We can compute them using image brightness E(x,y,t), and its spatial and time deriva- 
tives, E r and Et, over some region / in the image plane. Since there are only eight motion 
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and surface parameters to recover (There are three components of each of oj, t, and n. 
However, as mentioned earlier, the translational velocity and the surface normal can be 
recovered up to a scale factor.), only eight of the p,'s are independent. This implies that 
we can arbitrarily fix one of the essential parameters, and compute the remaining eight 
using eight independent points (At each point, we get one constraint and we have eight 
unknowns to recover). 

Let p' = (pi,P2>— >P8>0) denote the solution obtained by setting the last element 
equal to zero. If we define: 

P , = (Pi>P2'—'P8)> and a= (ri5i,ri52,riS3,r25i,r252,r 2 S3,r35i,r352) r , 

then the above constraint equation reduces to: 

a T p' + c = 0. 

Using eight non-colinear points (this is necessary to guarantee that the resulting equations 
are independent) , we can solve the following linear matrix equation: 

Ap' + c = 0, 

where: 

A = (ai,a2,...,ag) , c = (ci,C2,...,C8)- 

The solution of the above equation is: 

p' = -A _1 c. 

Image brightness values are distorted with sensor noise and quantization error. These 
inaccuracies are further accentuated by methods used for estimating the brightness gra- 
dient. Thus it is not advisable to base a method on measurements at just a few points. 
Instead we propose to minimize the error in the brightness constraint equation over the 
whole region I in the image plane. So we choose the essential parameters that minimize: 



//, 



;T~; 



(a p + c) dxdy. 
I 

The solution, in this case, is given by: 

p' = — ( / / aa T dxdy)~ 1 (l I ckdxdy). 

Note that, in general, the true pg is nonzero. We can show that the solution obtained 
through the assumption that p' 9 = 0, p', and the true solution (denoted by p) are related 
by the equation: 

p = p' + p 9 u, u = (1, 0, 0, 0, 1, 0, 0, 0, 1) T . 

The proof goes as follows. Since s = (E r x z) x r, then r T s = r • ((E r x z) x r) =0. For 
any arbitrary constant /, such that L = /I (I is the identity matrix), we have: 

r r Ls = 0. 
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If we add this to our costraint equation: 

c + r T (W + nt T + L)s = 0. 
It is immediately apparent that any P' of the form: 

( Pi P2 P3\ 
P'= P\ Ps Pel =W + nt T + L 

V?7 Ps PqJ 

will also satisfy our constraint equation. Therefore, the two solutions for the essential 
parameters are related by: 

P = P' - L, 

or in terms of p and p': 

p = p'-/u, u= (1,0,0,0,1,0,0,0,1), 
for some constant /.It follows that pg = p' g — I. Now if p' Q — 0, then / = — pg, and so: 

P = P' + P9U- 

In term of matrices P, and P': 

P = P' + p Q I. 

The procedure for determining pg, and subsequently P is presented in the appendix. For 
now, we assume that P is known (note that we have, so far, shown how to compute P' 
and not P). 

4. Recovering Motion and Structure 

We now show how to compute the parameters of the translational motion and the plane 
orientation from the essential parameters. Once these are known, the rotational param- 
eters are determined from: 

W = P - nt r 

Since W is skew-symmetric, and from the definition of P: 

P* = P + P T = {W + W T ) + (tn T + nt T ) = (tn T + nt T ). 

Consider the special case where t || n (n = A;t, for some constant A;). In this case, 

P* = 2ktt T 

Computing t from the above equation, and subsequently n is a simple practice of algebra. 
Therefore, from now, we assume that t x n / 0. 

Since the translational vector and the surface normal can be recovered up to a scale 
factor, we can, without loss of generality, let t r t = 1. We will present our results in 
terms of the normalized n and t, that is: 

„ _ n ~_ t 

n -R * - i*r - *• 
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Note that: 

TraceP* = (n T t + t T n) = 2(n • t) = 2 |n| (n • t), 

and therefore, we can always copmute n from: 

. . „ , . TraceP* 

n = |n|n ' |n| = ^TT 

We also have to normalize P*. From now, P* denotes the normalized matrix, so that: 

P* = (tn T + iit T ). 

In order to present the solutions for n and t, it is necessary to express the eigenvalue 
decomposition of P* in terms of these vectors. We will do so in the next lemma. 

Lemma 1: Let P* = XJAU T be the eigenvalue decomposition of P* = (tn r + nt r ), 
and let r = |TraceP*. Then: 

/r-1 \ 
A= , 

V r + 1/ 

t — n txfi t + ii 



U = 



W2{1-t)) |txn| \/2(l + r)) 



Proof: (tn T + nt T )txn = t(fi-(txn)) +ii(t-(txn)) = = 0-(txn), and so (txn) 

is the eigenvector corresponding to the zero eigenvalue of P*. Since it is symmetric, P* 
has 3 orthogonal eigenvectors. The other two eigenvectors are, therefore, in the plane 
containing t and ii. Let u = at + /3h and A denote an eigenvector-eigenvalue pair for 
some a and (3 (to be determined). Then 

(tn r + nt T )(at + (3h) = A(at + /m), 

which simplifies to 

[a(t • n) + /3{n ■ n)]t + [a(t • t) + /?(t • n)]n = Aat + A^n. 

Since (t • n) = r, we have: 

(r,j (,-)-■• 

The solutions for A are given by: 

A = r±l. 

Substituting for A into the earlier equations, we get: 

a = ±0. 
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Using these into the equation for u, and normalizing the results yield: 

tin 

u = . =. 

V2{1±t) 

Note that since \t\ < 1, Ai < A 2 = < A3. 1 

We can now detrmine t and n. Let u t - denote the i-th column of U. From lemma 1: 

t—h t+h 



%/2(l-r)' V2(1 + t) 

From the expressions for the eigenvectors, it follows that: 



t + n = \/2(l + r)u 3 , t - n = \/2(l - r)ui 
Solving these for t and n: 



1-T 
— 7T- u l> 



Since the choice of the signs of the eigenvectors are arbitrary, we should repeat the above 
procedure with the sign of either Ui or 113 reversed: 

t — n t + fi 

ui = u 3 



V2(l-r)' \/2(l + r)" 

In this case, the second solution is obtained from the following equations 



T 

-ui, 



Note that this is the dual solution for the planar surfaces. The special case of t || n 
corresponds to when the matrix P* has multiple eigenvalues. Then, either: 

1) r = 1, for which Ai = A2 = 0. Then the two solutions merge to the single one: 

n = t = u 3 , 

or: 

2) t = — 1, so that A2 = A3 = 0. Here, the unique solution is given by: 

n = — t = ui. 

As mentioned earlier, we can determine n through proper scaling of n, and solve for 
the rotation parameters by substituting the solutions for n and t into the equation: 

W = P - nt T . 
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Even though we gave a complete and compact proof of the dual solution in an earlier 
paper[15], it is intriguing to confirm those results with our closed form solution. We 
showed that the two solutions are related by: 

n* = A;t, t* = 1/k n, uj* = u + n x t, 

where k is arbitrarily chosen to scale the translation vector and plane normal. The two 
solutions given in lemma 2 for n and t already satisfy the duality relationship given 
above. It remains to show the same for the rotation parameters. We only have to show 
that: 

W* + n*t* T = P, 

where 

/ nit 2 — nzh n\tz — n 3 <i A 



i U2H — nii2 u 7*2*3 — '*3*2 j > 

V n 3 ti - nit 3 n 3 i 2 - n<>tz J 



W = W + | nitx - n\t% n 2 tz - 7*3*2 

l -2 — ^2*3 
T + „T 



n 3 ii — n\tz n 3 i 2 — ^2*3 
or 

W* = W + nt J - tn 1 . 



Substituting for n*, t*, and W* into the earlier equation, and simplifying the results, we 
get: 

W + nt T = P. 

5. Summary 

In this paper, we presented a closed form solution for recovering the motion of an 
observer with respect to a planar surface without having to compute the optical flow 
as an intermediate step. We only need the image intensity gradients at 8 non-colinear 
points, however, in general, it is better to compute these gradients in a larger portion 
of the image to reduce the noise effects. We first employed a constraint equation we 
developed for planar surfaces to compute 9 intermediate parameters, the elements of a 
3x3 matrix. We referred to them as essential parameters. The special structure of this 
matrix allows us to compute the motion and plane parameters very easily. 
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Appendix 

In the previous sections, we showed how the motion paprameters can be recovered 
once the essential paprameters are known. However, the brightness change constraint 
equation allowed us to determine the matrix P' (a particular solution of P with the last 
element set to zero). We showed that the two solution are related by: 

P = P' + p 9 I. 

Here, we show how to determine pg, and consequently P. For simplicity, let pg = I, and 
let P = UAV r denote the eigenvalue decomposition of P, where (A = diag(Ai, A2, A3)), 
then: 

P' = UAV T - II = UAV T - ZUV r . 

If L = II then: 

P' = UAV r - ULV r = U(A - L)V r . 

Similarly, if P* = P + P r = UAU r denotes the eigenvalue decomposition of P*, then: 

P'* = p' + p^ = p + p^ - 2L = P* - 2L = U(A - 2L)U T . 

In lemma 1, we showed that Ai<A2 = 0<A3 (and when t || n, we get two zero 
eigenvalues) . Therefore, the eigenvalues of P'* can be arranged in the form: 

Ai - 21 < -21 < A3 - 21. 

It follows that / = -iA 2 (P'). 

So in summary, we assume that pg = 0, and solve for the essential parameters (ele- 
ments of P'). The eigenvalue decomposition of P'* allows us to determine the unknown 
shift (pg = — IA2), and then, P from: 

P = P'-iA 2 (P')I 

Finally, the solution of the motion and structure parameters are determined from the 
equations given earlier in terms of the trace and eigenvectors of P* (note that the eigen- 
vectors of P* and P'* are the same). 
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